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1.  OBJECTIVES 

The  principal  goal  of  this  project  was  to  determine  the  usefulness  of  a  specific  tech¬ 
nique  for  deconvolution,  minimum  entropy  deconvolution  (MED),  in  detecting  seismic 
source  multiplicity  and  its  potential  for  discriminating  ripple-hred  explosions  from  other 
seismic  events.  The  technique’s  usefulness  in  detecting  multiple  sources  which  are  not  uni¬ 
formly  separated  in  time  and  are  not  easily  detected  by  spectral  methods,  has  been  investi¬ 
gated. 

This  report  consists  of  a  brief  introduction,  an  overview  of  the  MED  algorithm,  a 
description  of  tests  of  the  algorithm  on  synthetic  data,  the  results  of  applying  the  method 
to  data  from  known  and  suspected  quarry  blasts,  and  a  summary  of  the  method’s 
effectiveness. 

2.  INTRODUCTION 

Existing  techniques  for  discriminating  ripple-fired  quarry  explosions  from  other 
seismic  sources  rely  on  spectral  characteristics.  Smith  and  Grose  (1987)  observed  that 
under  certain  conditions  ripple-firing  can  produce  high-firequency  spectral  peaks  in  the  P 
spectra.  Several  investigators,  including  Baumgardt  and  Ziegler  (1988),  Stump  and  Reamer 
(1988),  Smith  (1989),  and  Hedlin  et  al.  (1989),  have  observed  spectral  modulation  in  the 
coda  of  data  from  mine  explosions.  These  spectral  characteristics  can  be  effectively 
modeled  by  assuming  that  waveforms  from  individual  shots  superpose  linearly  to  produce 
the  observed  signal.  Stump  and  Reinke  (1988)  present  experimental  evidence  to  support 
the  validity  of  linear  superposition.  Orcutt  and  Hedlin  (1991)  have  presented  an  automated 
algorithm  for  discriminating  ripple-fired  events,  which  relies  on  the  time-independence  of 
the  spectral  modulation  in  the  coda. 

The  spectral  modulations  produced  by  ripple  firing  are  greatly  reduced  when  the  time 
separations  between  the  individual  shots  of  a  ripple-fired  explosion  are  not  uniform.  The 
work  of  Stump  and  Reamer  (1988)  indicates  that  actual  shot  separation  times  can  be 
different  from  the  designed  separation  times  by  as  much  as  34%.  The  motivation  behind 
the  current  project  was  to  find  a  time-domain  discriminant  which  could  complement  spec¬ 
tral  discriminants  in  the  case  of  irregular  shot  times.  The  MED  algorithm  does  not  depend 
on  regular  shot  intervals,  and  theoretically  could  deconvolve  multiple  shot  times  when  the 
associated  spectral  modulations  are  very  weak.  In  practice  though,  there  are  several  factors 
which  limit  the  usefulness  of  MED  as  a  discriminant  for  ripple-fired  explosions.  These  are 
discussed  below  after  a  brief  overview  of  the  MED  algorithm. 

3.  MINIMUM  ENTROPY  DECONVOLUTION 

First  introduced  by  Wiggins  (1977,  1978),  minimum  entropy  deconvolution  is  a 
method  of  generating  a  linear  filter  capable  of  nansforming  its  input  signal  into  a  signal 
with  its  energy  concentrated  into  as  few  spikes  as  possible.  This  is  referred  to  as  maximiz¬ 
ing  the  spikiness  of  the  signal,  which  is  equivalent  to  minimizing  its  disorder,  or  minimiz¬ 
ing  the  entropy.  If  the  input  signal  can  be  characterized  as  the  convolution  of  a  seismic 
source  wavelet  with  a  series  of  spikes,  MED  represents  a  deconvolution  operation.  We 
assume  that  the  seismic  signal  from  a  multiple  shot  blast  can  be  modeled  as  the 
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superposition  of  time  delayed  replications  of  the  signal  from  a  single  shot.  This  implies 
that  a  deconvolution  method  could  be  used  to  extract  the  shot  delay  times  from  the  multi¬ 
ple  shot  signal.  The  design  of  the  MED  filter  centers  around  a  mathematical  measure  of 
the  spikiness  of  a  signal.  For  example,  the  varimax  norm  suggested  by  Wiggins  (1978)  is 


I.yf 


v  = 


,2 


(Zyf) 


(1) 


where  y,-  represents  the  seismic  signal.  Maximizing  V  with  respect  to  the  filter  coefficients 
leads  to  an  equation  which  can  be  iteratively  solved  for  the  MED  filter  coefficients.  The 
solution  of  the  equation  is  particularly  simple,  since  it  involves  a  Toeplitz  autocorrelation 
matrix.  Filter  designs  based  on  other  norms  have  also  been  suggested.  Ooe  and  Ulrych 
(1979)  investigated  the  use  of  an  MED  norm  including  an  exponential  transformation  with 
a  damping  coefficient  to  control  the  effect  of  noise.  The  transformation  is  defined  by 


yf 


zi  =  \  -  exp(  — -  ) 


2s^ 


(2) 


where  s  is  a  constant.  The  Ooe-Ulrych  norm  is  defined  by 


i 


(3) 


We  have  experimented  with  both  of  these  norms  and  have  found  the  latter  one  (3)  to 
be  the  more  useful.  The  iterative  solution  for  the  filter  coefficients  which  maximize  either 
norm  (1)  or  norm  (3)  is  a  linearization  of  a  nonlinear  equation,  and  thus  the  solution  can 
vary  greatly  depending  on  the  initial  conditions  and  the  value  of  the  constant  s.  In  addition 
to  the  complication  of  nonuniqueness,  we  have  found  from  tests  with  synthetic  data,  that 
the  filter  which  actually  maximizes  a  particular  norm  frequently  does  not  deconvolve  the 
data  as  well  as  one  of  the  filters  from  an  earlier  step  in  the  iterative  solution.  Therefore,  in 
practice  one  is  forced  to  review  the  effectiveness  of  each  filter  in  the  iteration.  The  param¬ 
eter  s  controls  the  sensitivity  to  small  features  in  the  data.  As  we  shall  demonstrate  below, 
in  this  application  of  the  MED  algorithm,  the  results  are  very  sensitive  to  the  value  of  s. 
Additional  parameters  in  the  problem  are  the  length  of  the  filter  and  the  position  and 
length  of  the  input  data  window. 

4.  TESTS  WITH  SYNTHETIC  DATA 

The  MED  algorithm  was  developed  as  a  tool  for  processing  reflection  seismic  data. 
The  objective  in  that  case  is  to  deconvolve  the  impulse  response  of  the  transmission  path. 
In  general,  the  time  separation  between  impulses  (the  travel  time  between  reflective  boun¬ 
daries)  is  greater  than  the  dominant  period  of  the  data.  In  our  application  of  the  MED 
technique,  the  time  separation  between  impulses  (ripple-fired  shots)  can  be  much  smaller 
than  the  period  of  the  data.  Therefore,  in  our  application,  the  deconvolution  relies  more 
heavily  on  high  frequency  components  of  the  data.  In  this  section  we  discuss  several 
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examples  of  applying  the  algorithm  to  synthetic  data  to  illustrate  its  dependence  on  data 
sample  rate  and  frequency  content. 

Assuming  that  ripple-fired  shot  delay  times  are  approximately  25msec,  we  have  deter¬ 
mined  the  data  requirements  for  a  successful  deconvolution.  We  begin  first  with  an  ideal¬ 
istic  example  of  synthetics  with  a  very  high  sample  rate,  in  which  the  MED  algorithm 
works  well  and  would  complement  spectral  methods  in  the  case  of  nonuniform  shot  delay 
times.  Then  we  illustrate  how  the  deconvolution  deteriorates  as  the  sample  rate  decreases 
or  the  period  of  the  source  wavelet  increases. 

Figure  la  shows  a  20Hz  wavelet  sampled  at  4(X)sps,  which  was  convolved  with  the 
uniform  impulse  series  (25msec  delay  time)  of  Figure  lb  to  produce  the  synthetic 
waveform  shown  as  Figure  Ic.  The  same  wavelet  was  convolved  with  the  nonuniform 
impulse  series  of  Figure  Id  to  generate  the  waveform  shown  as  Figure  le.  The  amplitude 
of  the  impulses  in  Figure  Id  has  a  10%  random  variation  and  the  delay  times  between 
impulses  are  randomly  distributed  between  15  and  35msec  (25msec  +/-  10msec).  The 
spectra  of  both  synthetic  waveforms  is  shown  in  Figure  2,  to  illustrate  the  effect  of  nonun¬ 
iform  shot  separation  times.  The  spectrum  (Figure  2a)  of  the  waveform  produced  from 
the  uniform  impulse  series  has  prominent  peaks  at  multiples  of  40Hz,  corresponding  to  the 
uniform  25msec  impulse  separation.  In  contrast,  the  spectrum  (Figure  2b)  of  the  waveform 
produced  from  the  nonuniform  impulse  series  does  not  have  any  modulations. 


Figure  1.  a)  A  20  Hz  wavelet  sampled  at  400sps.  b)  A  series  of  39  uniformly 
spaced  impulses  separated  by  25msec.  c)  The  wavelet  "a"  convolved  with  the 
impulse  series  "b".  d)  A  series  of  39  impulses  with  separations  randomly  distri¬ 
buted  between  15  and  35msec.  e)  The  wavelet  "a"  convolved  with  the  nonuniform 
impulse  series  "d". 
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Figure  2.  a)  The  spectrum  of  the  waveform  show  in  Figure  Ic,  clearly  showing 
modulations  due  to  the  uniform  impulse  separation,  b)  The  spectrum  of  the 
waveform  show  in  Figure  Id,  for  the  case  of  nonuniformly  spaced  impulses. 


Both  of  the  waveforms  shown  in  Figures  Ic  and  Id  can  be  deconvolved  with  the 
MED  algorithm  after  highpass  filtering  at  50Hz.  The  deconvolution  of  the  nonuniform 
impulse  series  is  shown  in  Figure  3.  Figure  3a  shows  the  synthetic  waveform  from  Figure 
le  highpass  filtered.  The  deconvolved  waveform  is  shown  as  Figure  3b,  with  the  input 
nonuniform  impulse  series  below  it  for  comparison.  The  deconvolution  (Figure  3b)  recov¬ 
ers  the  impulse  times  very  well.  It  does  not  recover  the  amplitudes  of  the  input  impulses 
as  well,  but  the  deconvolution  would  certainly  be  a  clear  indication  of  source  multiplicity. 

If  the  period  of  the  source  wavelet  (Figure  la)  is  increased,  or  the  sample  rate  of  the 
synthetic  waveform  (Figure  le)  is  decreased,  the  effectiveness  of  the  MED  algorithm  is 
diminished.  Figure  4  illustrates  the  effect  of  increasing  the  period  of  the  source  wavelet. 
Three  deconvolutions  are  shown  for  source  wavelets  with  periods  of  approximately 
SOmsec  (20H2),  67msec  (15Hz),  and  100msec  (lOHz).  The  MED  algorithm  is  still  able  to 
recover  all  of  the  impulses  for  the  15Hz  wavelet,  whose  period  is  about  2.5  times  the 
average  impulse  separation  of  25msec.  The  deconvolution  with  the  lOHz  wavelet  fails  to 
recover  all  of  the  impulses  and  produces  a  few  spurious  spikes. 
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Figure  3.  a)  The  synthetic  waveform  of  Figure  le  highpass  filtered  at  50  Hz.  b) 
The  deconvolution  of  the  waveform  "a",  c)  The  nonuniform  impulse  series  which 
was  used  to  generate  the  synthetic,  shown  for  comparison. 


0  1 

seconds 


Figure  4.  Deconvolutions  of  three  synthetics  made  with  three  different  source 
wavelets  of  varying  period  and  the  nonuniform  impulse  series  (Figure  Id),  a) 
20Hz  wavelet,  b)  I5Hz  wavelet,  c)  lOHz  wavelet,  d)  Deconvolution  using  20Hz 
wavelet,  e)  Deconvolution  using  15Hz  wavelet,  f)  Deconvolution  using  lOHz 
wavelet. 
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Figure  5  illustrates  the  effect  of  decreasing  the  sample  rate  of  the  synthetic  waveform 
before  applying  the  deconvolution  technique.  Four  deconvolutions  of  the  synthetic  shown 
in  Figure  le  (generated  with  the  20Hz  wavelet  and  the  nonuniform  impulse  series  with  a 
2Smsec  average  impulse  separation)  are  shown  for  sample  rates  of  400,  325,  250,  and  100 
samples  per  second.  The  best  result  is,  of  course,  obtained  with  the  400sps  synthetic.  The 
algorithm  works  almost  as  well  at  325sps.  It  begins  to  fail  at  250sps  and  loses  most  of  its 
resolution  at  150sps,  recovering  only  a  few  impulses. 


seconds 


Figure  5.  Deconvolution  of  synthetics  made  with  the  20Hz  wavelet  (Figure  la)  and 
the  nonuniform  impulse  series  (Figure  Id),  using  different  sample  rates,  a)  400sps. 
b)  325sps.  c)  250sps.  d)  ISOsps. 


A  final  test  of  the  MED  algorithm  on  synthetic  waveforms  was  made  to  determine 
the  effectiveness  of  the  deconvolution  of  long  data  windows.  Synthetic  waveforms  were 
produced  by  summing  time  delayed  copies  of  recorded  data.  An  example  of  a  successful 
deconvolution  of  this  type  of  synthetic  is  shown  in  Figure  6.  Figure  6a  displays  20 
seconds  of  a  short  period  record  (60  sps),  which  was  convolved  with  an  impulse  series 
(Figure  6b)  consisting  of  27  impulses  with  a  random  separation  between  40  and  50msec, 
to  produce  the  synthetic  waveform  shown  as  Figure  6c.  The  MED  algorithm  was  able  to 
recover  about  a  third  of  the  impulses  (Figure  6d).  The  algorithm  was  not  able  to  decon¬ 
volve  similar  synthetics  made  with  60sps  data  and  time  delays  less  than  50msec. 
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Figure  6.  Deconvolution  of  a  synthetic  made  by  summing  recorded  data,  a)  A  short 
period  record  (60sps).  b)  Impulse  series  consisting  of  27  impulses  with  random 
separation  between  40  and  50msec.  c)  The  record  "a"  convolved  with  the  impulse 
series  "b".  d)  deconvolution  showing  about  9  recovered  impulses  and  two  spurious 
spikes. 

The  examples  discussed  above  place  some  bounds  on  the  effectiveness  of  the  MED 
algorithm  in  terms  of  data  sample  rate,  data  frequency  content  and  data  window  length,  for 
a  simplistic  representation  of  the  seismogram.  As  might  be  expected,  we  also  found  the 
technique  to  be  ineffective  in  the  presence  of  a  small  amount  high  frequency  noise.  In 
addition,  the  algorithm  is  quite  sensitive  to  all  of  its  parameters:  filter  length,  damping 
coefficient,  data  window  length  and  the  number  of  iterations.  The  successful  deconvolu¬ 
tions  shown  were  achieved  only  by  time  consuming  interactive  experimentation  with  these 
parameters. 

5.  APPLICATION  TO  DATA 

We  have  tested  the  MED  algorithm  on  data  from  known  and  suspected  quarry  blasts. 
Figure  7  shows  a  vertical  component  short-period  record  from  station  KK  of  the 
Soviet/NRDC  data  set,  at  a  distance  of  180km  from  the  source.  Many  events  from  this 
data  set  are  suspected  to  be  mining  blasts,  and  its  high  sample  rate  (250sps)  made  it  an 
appropriate  choice  for  attempting  to  deconvolve  sources  separated  by  time  delays  greater 
than  approximately  25msec.  The  deconvolution  results  shown  in  Figure  7  are  the  most 
interesting  results  obtained  after  examining  more  than  50  records  from  this  data  set.  The 
middle  waveform  in  Figure  7  is  the  first  6  seconds  of  the  signal  shown  on  the  top.  A  2  8- 
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second  (7(X)  point)  MED  filter  was  generated  for  this  data  window.  The  filtered  or  decon¬ 
volved  data  is  shown  on  the  bottom  of  Figure  7.  The  output  of  the  filter  operation  is  typi¬ 
cally  time  shifted  by  an  arbitrary  amount,  because  the  MED  algorithm  does  not  constrain 
the  phase  of  the  filter.  The  deconvolved  data  consists  of  three  pairs  of  impulses  separated 
by  380  and  170  msec.  The  time  between  the  impulses  of  each  doublet  is  a  nearly  constant 
47  msec.  We  suspect  that  the  three  doublets  represent  three  different  phase  arrivals.  The 
doublet  nature  could  be  due  to  source  multiplicity,  although  this  was  the  only  record  from 
this  data  set  for  which  we  obtained  such  a  result. 


0  20  40  60 


seconds  from  1987079  08:08:7 


Figure  7.  top:  Vertical  component  record  from  station  KK  for  an  event  at  (49.9N, 
73. IE),  middle:  Data  window  at  the  beginning  of  the  signal  for  which  an  MED 
filter  was  generated,  bottom:  The  deconvolved  data  window.  The  time  between  the 
impulses  in  each  of  the  three  doublets  is  a  nearly  constant  47  msec. 
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Figure  8.  top:  Vertical  component  borehole  record  of  a  quarry  blast  at  a  distance  of 
32  km.  middle:  Data  window  at  the  beginning  of  the  signal  for  which  an  MED 
filter  was  generated,  bottom:  The  deconvolved  data  window. 


The  algorithm  has  also  been  tested  on  data  from  a  known  quarry  blast,  for  the  case  of 
a  shot  separation  time  smaller  than  the  data  sampling  interval.  For  a  quarry  blast  about  32 
km  from  the  Oklahoma  Geophysical  Observatory,  the  average  delay  time  between  shots 
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was  approximately  9  msec,  while  the  sample  rate  for  the  record  from  LNO  is  60  sps  (17 
msec  sample  interval).  There  were  a  total  of  27  shots  spanning  approximately  .25  seconds, 
arranged  in  two  rows.  The  delay  time  from  the  start  of  the  first  row  to  the  start  of  the 
second  row  was  42  ms.  Figure  8  shows  a  vertical  component  borehole  recording  of  the 
signal  from  the  quarry  blast  at  station  LNO.  The  middle  trace  of  Figure  8  shows  the  first  4 
seconds  of  the  signal  that  was  input  to  the  MED  algorithm.  A  4-second  MED  filter  was 
generated  which  deconvolved  this  signal  into  the  trace  shown  at  the  bottom  of  Figure  8.  It 
is  tempting  to  interpret  the  deconvolved  trace  as  4  phase  arrivals,  the  first  three  consisting 
of  two  closely  spaced  spikes  due  to  source  effects.  The  time  separation  between  the  spikes 
in  each  doublet  is  50  msec.  Could  the  double  nature  of  the  deconvolution  be  caused  by  the 
delay  time  between  start  of  each  row  of  shots? 

Because  of  the  nonlinear  behavior  of  the  algorithm,  it  is  very  difficult  to  have 
confidence  in  such  an  interpretation  of  the  deconvolution.  Different  deconvolutions  can  be 
obtained  for  the  same  data  window  by  simply  varying  the  filter  length  or  damping 
coefficient.  Figure  9  shows  four  deconvolutions  for  the  data  window  of  Figure  8.  Each 
trace  in  Figure  9  was  generated  by  simply  varying  the  filter  length  or  the  damping 
coefficient.  All  four  deconvolutions  have  some  doublet  impulses  and  some  of  the  impulses 
are  common  to  all  of  the  traces.  Without  prior  knowledge  of  an  expected  result,  as  was  the 
case  with  the  synthetic  waveforms,  the  nonuniqueness  of  the  MED  deconvolution  clearly 
reduces  its  usefulness. 


0  12  3  4 


seconds 

Figure  9.  Four  deconvolutions  obtained  for  the  data  window  shown  in  the  Figure  8. 
Each  trace  corresponds  to  a  different  value  for  the  filter  length  and/or  the  damping 
coefficient. 
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CONCLUSIONS 

Our  experience  with  applying  this  deconvolution  method  to  synthetic  and  real  data 
from  quarry  blasts  leads  us  to  several  observations  about  the  nature  of  the  technique.  In 
general,  the  algorithm  is  not  robust,  but  is  very  sensitive  to  the  filter  parameters,  filter 
length  and  filter  damping  coefficient.  It  is  also  sensitive  to  the  data  window  length  and 
position.  Tests  with  synthetic  data  indicate  that  the  filter  which  performs  the  best  deconvo¬ 
lution  frequently  does  not  correspond  to  that  which  maximizes  the  data  norm,  but  occurs 
earlier  in  the  iterative  solution  for  the  maximum.  The  method  is  dependent  on  high  fre¬ 
quency  data.  For  a  shot  interval  of  25msec,  our  experiments  with  synthetics  indicate  that 
data  above  50Hz  is  required  for  a  deconvolution.  For  cases  where  the  data  sampling  inter¬ 
val  is  adequate,  the  sensitivity  of  the  algorithm  will  remain  as  the  chief  obstacle  hindering 
attempts  to  develop  the  method  into  an  automatic  di.scriminant.  The  MED  algorithm  could 
be  a  useful  tool  for  deconvolving  overlapping  phase  arrivals,  for  interpretation  and  more 
precise  inter-arrival  timing. 
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